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Abstract 


We present a new single level effectively explicit implicit algorithm for 
gasdynamics. The method meets all the requirements for unconditionally 
stable global iteration over flows with mixed supersonic and subsonic zones in- 
cluding blunt body flow and boundary layer flows with strong interaction and 
sireamwise separation. For hyperbolic (supersonic flow) regions the method 
is automatically equivalent to contemporary space marching methods. For 
elliptic (subsonic flow) regions, rapid convergence is facilitated by alternat- 
ing direction solution sweeps which bring both sets of eigenvectors and the 
influence of both boundaries of a coordinate line equally into play. Point by 
point updating of the data with local iteration on the solution procedure at 
each spatial step as the sweeps progress not only renders the method single 
level in storage but, also, improves nonlinear accuracy to accelerate conver- 
gence by an order of magnitude over related two level linearized implicit 
methods. The method derives robust stability from the combination of an 
eigenvector split upwind difference method (CSCM) with diagonally dominant 
ADI (DDADI) approximate factorization and computed characteristic bound- 
ary approximations. The properties and performance of the technique are 
demonstrated in a variety of quasi 1-D nozzle flows including completely 
subsonic or supersonic or mixed subsonic/supersonic with sonic points and 
shocks. The performance of the new method in 2-D and axisymmetric flows 
are compared against the precursor two level linearized scheme in an invis- 
cid inlet problem for which the exact solution is known and in a viscous 
transonic nozzle problem previously studied both experimentally and com- 
putationally. 
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I. Background 


The motivation for the present work stems from the following reason- 
ing. Parabolized Navier Stokes (PNS) procedures can effectively solve steady 
supersonic and boundary layer flow problems with small adverse pressure 
gradient at a fraction of the cost required by conventional time dependent 
procedures. But PNS procedures are ill-posed for subsonic and stream- 
wise reversed flow problems from the theoretical point of view, a failing 


that is manifested in grid sensitive numerical stability problems. On the 


other hand, certain eigenvector split upwind differenced implicit time depen- 
dent. procedures are unconditionally stable, practically robust and their ap- 
plication is not restricted with respect to the nature or type of compres- 
sible flow. A solution procedure that has the speed and the storage re- 
quirements of PNS schemes and the stability and applicability of the im- 
plicit upwind time marching schemes will be very useful and desirable. We 
give such a procedure here. Before we go on to describe this novel al- 
gorithm, we present a brief background to PNS and other related tech- 
niques. 

The Beam-Warming factored implicit algorithm 1 with the Baldwin- . 
Lomas thin layer viscous approximation 2 has provided the basis for two 
similar space marching (PNS) procedures 3,4 for the compressible Navier- 
Stokes equations. These PNS methods which are highly efficient from 
the points of view of both data storage and computer time have proven 
effective for flows 5,8 with favorable streamwise pressure gradient or with 
relatively small adverse pressure gradients. However, in the presence of 
strong adverse pressure gradient such as occurs in a wing or fin root regions 
the contemporary PNS methods suffer numerical stability problems and may 
infer streamwise separation even where separation doesn’t occur 7 . In such 
unseparated (perhaps weakly separated) regions, numerical stability may be 
maintained at the price of employing large amounts of artificial viscosity 
with a resulting loss in predictive accuracy and knowledge of the actual 
state of the flow. Where strong streamwise separation occurs the methods 
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are unstable and cannot proceed. Particularly for cue increasingly relevant 
'laminar flow situation that will be encountered at very high altitude by 
aerodynamic systems such as orbital transfer vehicles (OTV’s) and space 
shuttle, streamwise separation becomes a likely occurrence 8 in compression 
corners associated with canopies, pods, flared bodies, wing or fin roots and 
deflected control surfaces. Thus a more general technique is needed that 
is inherently stable for all types of upstream influence. At a minimum 
the mixed elliptic hyperbolic problem requires global iteration, preferably 
with type dependent differencing. Various steps in this direction have 
recently been taken by Rakich®, by Rubin and Reddy 10 and by Rizk and 
Chaussee 11 . 

Rakich® utilized global iteration with Vigneron’s Mach number depen- 
dent upwind and downwind splitting 4 of the pressure gradient for 
the streamwise momentum equation. The approach provides an im- 
provement in both accuracy and stability for both weakly interactive 
boundary layer flow and strongly interactive flow without (significant) 
streamwise separation. Global iteration is carried out by march- 
ing the PNS equation repeatedly only along the downstream direc- 
tion. 


For incompressible flow Rubin and Reddy 10 (also Lin and Rubin 12 
for supersonic flow) introduced type dependent (upwind) differencing 
of the streamwise velocity along with pressure splitting in an im- 
plicit method involving a staggard grid - dependent variable location 
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scheme. This method was extended to compressible flow by Khosla 
and Lai 13 . The method admits strong interaction with streamwise 
separation but is not homogeneous across the spectrum of Mach num- 
ber regimes of interest here. Other PNS related methods for subsonic 
flow and boundary layer flow with streamwise reversal have recently 
been reviewed by Brown 14 in conjunction with a new staggered grid 
scheme. 

Most immediately relevant, Rizk and Chaussee 11 presented a hybrid 
technique of space marching in supersonic zones and relaxation with a 
Beam-Warming time dependent central difference method over zones of 
strong upstream influence with streamwise separation. They have two 
variations of the time dependent method: one fully implicit requiring two 
levels of storage and matrix inversion procedures in all space directions, 
the other explicit in all space directions but the thin layer direction and 
requiring only one level of storage as in a space marching algorithm. 
Both procedures require the same several hundred iterations minimum to 
reasonably converge, with the semi explicit procedure requiring substantially 
less machine time. However, both procedures are inherently unstable, 
except for the use of artificial dissipation. The marginal stability coupled 
with the lack of type dependent differencing and well posed boundary 
approximations all contribute to slow convergence. In balance, while 
workable, the hybrid approach leaves something to be desired from the points 
of view of convenience, computer time and internal consistency of the global 
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solution procedure. 

\ 

A new globally iterated scheme related to that which will be presented 
here, has been presented by Moretti 15 . In his procedure for steady invis- 
cid flows, the Euler equations are cast in Riemann variables and the result- 
ing uncoupled equations are solved by integrating each Riemann variable 
separately, sweeping back and forth alternatively along each coordinate line. 
The coupling between the equations and, thus, the non-linearity are intro- 
duced only through the boundaries and the updating of state after com- 
plete sweeps. As can be inferred from the results to be presented here, the 
reduced coupling inherently limits the rate of convergence of Moretti’s method 
and the procedure is- not extendable to the compressible Navier-Stokes equa- 
tions. 

II. New Universal Single Level Scheme CSCM-S 

The CSCM flux difference eigenvector split upwind implicit method 18,17,18 
for the inviscid terms of the compressible Navier-Stokes equations provides 
the natural basis for an unconditionally stable space marching tech- 
nique through regions of subsonic and streamwise separated flow. In 
such regions the split method can be likened to stable marching of 
each scalar characteristic wave system in the direction of its associated 
eigenvalue (simple wave velocity). In supersonic flow, where all eigen- 
values have the same sign, the method automatically becomes similar 
to the referenced PNS techniques based on the Beam-Warming fac- 
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tored implicit method with the Baidwin-Lomax thin layer viscous ap- 
proximation. 

% 

Compared to contemporary central difference methods, the CSCM 
characteristics based upwind difference approximation with its inherent 
numerical stability leads to greatly reduced oscillation and greater ac- 
curacy in the presence of captured discontinuities such as shocks, con- 
tacts and physical or computational boundaries. The correct mathemati- 
cal domains of dependence that correspond with physical directions of 
wave propagation are coupled with well posed characteristic boundary 
approximations 17 naturally consistent with the interior point scheme. The 
result is faster sorting out of transient disturbances and substantially 
more rapid convergence to the steady state. The splitting and the as- 
sociated time dependent implicit method have been described in detail in 
references (16) and (18) for quasi 1-D and 2-D planar or axisymmetric 
flow. 

In the following, we will sketch the differences between the time depen- 
dent method and the new space marching technique which we designate 
CSCM-S. The discussion will begin with the quasi 1-D inviscid formula- 
tion, present some results elucidating the properties and performance of the 
method, then give additional details entering into multidimensional inviscid 
and thin layer viscous procedures and, lastly, present early 2-D solutions ob- 
tained with the new single level scheme in problems solved previously 18 with 
the time dependent method. 
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Quasi 1-D Formulation 


\ 

The general jth interior point difference equations for the time dependent 
v. CSCM upwind implicit method is -written 

(/ + A+ V + A~A)6qj = -A+ Aq) - A~ Aq) n . (1) 

J ^ J 

where V and A are backward and forward spatial difference operators. In 
the notation the interval averaged matrices between node points j and j + 1 
are tabled j. The right hand side of equation (1) is written for the first order 
method. Higher order methods in space are given with results in references 
16 and 18. 

Central to its accurate shock capturing capability, the CSCM con- 
servative flux difference splitting has the “property U” put forth by 
Roe 10 

( 2 + +2~)A?) j 3AF) j = F i+l - Fj (2) 

Here q is the conservative dependent variable vector and F is the as- 
sociated flux vector. The matrices A + and A - are the splittings of 
the CSCM interval averaged Jacobian matrix according to the signs of 
the averaged eigenvalues. Thus in the equation for the jth grid point, 
A + A^g)y_i represents stable characteristic spatial differencing backward 
for positive eigenvalue contributions and A~ A^q)j, forward for negative 


ones. 
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ORIGINAL PAGE 13 
OF POOR QUALITY 

With 5? = ? n + 1 — g n , equation (1) defines a two level linearized coupled 
t)lock matrix implicit scheme that can be solved by a block tridiagonal 
procedure. In reference (18) a new (DDADI) approximately factored alter- 
nating sweep bidiagonal solution procedure for equation (1) is presented that 
is shown to be very robust and is effectively explicit, i.e. requires only a 
decoupled sequence of local block matrix inversions rather than the solution 
of the coupled set. For the forward sweep the bidiagonal solution procedure 
can be written 


(I+A+ -A-)5q\=RHS-\-A + 6q\_ l (3) 

For the linear problem, i.e. constant coefficient case of stability 
analysis, equation (3) is equivalent to the single level space marching 
procedure 

(/+ 2+ -X+,” -2- a,)" (4) 

3 3 

Nonlinearity enters in the single level space marching form (4) in that 
at each step of the forward sweep the matrices A + are averaged 
between and q* rather than homogeneously at the old iteration 

level n. Similarly, companion backward space marching sweep that is 
symmetric to equation (4) and that is intimately related to the back- 
ward sweep of the alternating bidiagonal algorithm of reference (1<S) 
is 

U + — A~ )8qj =— £ + A -\-A ~ — 7C~ (5) 
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ORIGINAL PAGE fS 
OF POOR QUALITY 

As ■will be shown in a forthcoming paper 20 , the method given by equations 

(4) and (5) is von Neumann unconditionally stable for the scalar wave 
\ 

equation. The analysis shows the significance of DDADI approximate 
„ factorization in rendering both the forward and backward sweeps separately 
stable regardless of eigenvalue sign. Consequently as the local Courant 
number becomes very large, the robust method becomes a very effective 
(symmetric Gauss-Seidel) relaxation scheme for the steady equations, a fact 
which substantially contributes to the very fast performance that will be 
demonstrated. 

At a right computational boundary on the forward sweep we solve tne 
characteristic boundary point approximation 

(3f + -f )5q N — A + 9 jv ^ 

I 

q^ 1 = q{j and at a left, on the backward sweep 

(S 1 - - X~ to, =X-i n l -X- ?" + 1 (7) 

Following the solution of equations (6) and (7) the conservative state 
vector is iteratively corrected 17 to maintain the accuiacy of prescribed 
boundary conditions while not disrupting the representation of the com- 
puted characteristic variables running to the boundary from the inter- 
ior. Analysis of a model system with upwind differenced scalar equations 
and coupled boundary conditions was related to the linearized bidiagonal 
scheme 18 by Oliger and Lombard 21 ; the analysis also strongly supports the 
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numerically confirmed robust stability of the present nonlinear method for 
gasdynamics. 

With the updating at each step, "where in equation (5) 6qj = , 

it is clear that the symmetric pair of equations (4) and (5) serve to advance 
the solution two pseudo time (iteration) levels; whereas, the linear alternat- 
ing bidiagonal sweep algorithm of reference (18) advances the solution only 
one level. To maintain conservation to a very high degree, in single sweep 
marching in supersonic zones we iterate (at least) once locally at each space 
marching step. The local iteration serves to make the eigenvectors in the 
coefficient matrices consistent with the advanced state and, thus, provides 
improved accuracy for the nonlinear system. It appears effective to do this 
inner iteration everywhere, i.e. in both subsonic and supersonic regions, as 
the number of global iteration steps to convergence with two inner iterations 
has been found reduced by a factor of three to four. Since the computa- 
tional work per two steps is about the same for the single level and two level 
schemes and beyond the fact that one saves a level of storage in the space 
marching algorithm, the question arises: Can one get solutions in less com- 
putational work through faster convergence with the nonlinear space march- 
ing algorithm? 


One Dimensional Results 

First, we present results for supersonic flow with no shock in Shubin’s 
diverging nozzle. In purely supersonic zones, the experience with the 


10 


present method is that the solution can be marched accurately in one 
global iteration, as ought to be the case. Figure 1 shows the ex- 
act solution (in solid line) and the computed result from the first for- 
ward sweep. It is evident that the method correctly predicts the solu- 
tion to plotting accuracy in one forward sweep. With subsequent sweeps 
the error (the difference between the exact and the computed solution) 
reduces to machine accuracy in less than three global iterations. In 
fact, by increasing (from two) the number of inner iterations on the 
solution procedure at each space marching step one can guarantee con- 
vergence to prescribed accuracy in one forward sweep. This is also 
true cf contemporary locally linearized unsplit methods in supersonic 
flow. 

With the globally iterative nonlinear space marching formulation, 
early experience in two quasi 1-D nozzle problems with mixed supersonic- 
subsonic zones is that solutions are obtained in roughly an order 
of magnitude fewer iteration steps than had been required with the 
previously fast pseudo time dependent technique and block tridiagonal solv- 
ing. 


The two nozzle problems which are described and solved by Yee, 
Beam and Warming 22 and solved with the CSCM time dependent tech- 
nique in references (16) and (17) are Shubin’s diverging nozzle flow 
and Blottner’s converging-diverging nozzle flow. Both problems in- 
volve unmatched overpressures at the outflow which result in inter- 
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nal shock terminated supersonic zones and subsonic outflow. For 

the experiments involving flow of mixed type the same initial data 
given by Yee, Beam and Warming - a linear interpolation between 
inflow and outflow values for effectively exact solutions of the problems 
- is used that was used previously with the time dependent ap- 
proach. 

For flows of mixed type, in Figures 2 and 3 respectively, results 
are shown for successive forward and backward sweeps for five global 
iteration steps with Shubin’s and Blottner’s nozzle flows. In both cases, 
the exact solution as given by Yee, Beam and Warming is shown in solid 
line and the present computational results solved on a 51 point mesh, in 
boxes. Blottner’s nozzle flow is shown converged after 10 global iteration 
steps. There is substantial evidence in other results not shown that 
with further work the number of global iterations required to compute 
flows such as Blottner’s will be reduced by a factor of two, to about 
five. 


In Figure 4, we show a subscritical, i.e. complete!}' subsonic, flow 
solution computed in only two global iteration steps for the Blottner nozzle 
geometry with different inflow conditions. Here the exact analytical solution 
derived by Venkatapathy is shown in solid line and our computed results in 
boxes. 


The alternating direction sweeps in our method have been derived 
directly out of theory for solving the implicit set of difference equations. 



However, with a moments thought, one can see mechanistically, numeri- 
cally speaking, that omitting the backward sweep from the pair and 
globally iterating only with the forward sweep equation (4) will result in 
permitting the influence of a subsonic outflow boundary (or interior dis- 
turbance) to propagate upstream only one grid point per global itera- 
tion. In such a case, which relates to other global iteration methods 
found in the literature and that sweep only in the main flow direction, 
rate of convergence is greatly inhibited relative to symmetric sweeping 
by a factor of order roughly the number of grid points in the subsonic 
zone. Mathematically, this inhibition is the result of the failure to in- 
clude the effect of the eigenvectors governing upstream influence in the 
implicit process but to treat these waves explicitly with effective CFL 
unity. 

In Figure 5 we illustrate the progress of the transient solution to 
the subcritic al nozzle problem after 15 forward sweeps, with the back- 
ward sweeps omitted. One can clearly see that the wave influence of 
the outflow boundary has progressed only 15 mesh points forward of the 
outflow boundary. In Figure 6 the transient solution is shown after 60 
steps which is beyond one characteristic transit time (equivalent to 50 
mesh intervals) for the upwind wave to reach the inflow boundary. In 
Figure 7 we show the history of the RMS error in the primitive vari- 
ables. The solution is found to converge to roughly the same RMS er- 
ror after three characteristic times (150 steps) as the solution obtained 
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with the symmetric alternating sweep sequence after only 3 global itera- 
tions. 

Blottner’s supercritical nozzle problem which involves subsonic inflow ac- 
celerating through a sonic point to a supersonic zone terminated by a shock to 
subsonic outflow is the most computationally demanding of the test cases and 
indicates the capability for the method to compute simply and consistently 
over the subsonic forebody and base regions of blunt bodies in supersonic 
flow. Thus the need for separate time dependent codes will be obviated by 
this new method. 

Finally, in Figure 8 and 9, we present the convergence history for 
the present nonlinear scheme and the linearized time dependent scheme 
for completely subsonic and supersonic nozzle flows. The x-axis shows 
the number of iterations each scheme requires to reduce the exact error 
to five orders of magnitude for various Courant numbers. It is evident 
that the present scheme converges extremely fast at all CFL numbers 
compared with the method based on the linearized block tridiagonal sol- 
ver. 

Two Dimensional Formulation 

For two dimensional flow, assuming a marching coordinate £, inviscid 
terms 


y + v n +%-Ar, 


(6a) 
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and 

-%+A n q) k - l -y-A r ,q)k ( 66 ) 

are added to the left and right hand sides respectively of both the for- 
ward and backward sweep equations (4) and (5). For viscous flow, second 
centrally differenced, thin layer viscous terms are also added in the i] direc- 
tion as is conventionally practiced, e.g. Steger 23 . With the terms for 
the t) cross marching coordinate direction, the technique now becomes an 
implicit method of lines. Along each rj coordinate line, one can solve 
the equations coupled with a block tridiagonal procedure. Alternatively, 
a further DDADI bidiagonal approximate factorization can be employed 
in the 77 direction and solved either linearly as in reference (18) or non- 
linearly as here in the £ direction. As shown in the quasi 2-D numeri- 
cal experiments of reference (18). DDADI bidiagonal approximate factoriza- 
tion is stable for viscous as well as inviscid terms. Finally in reference 
(18) there is a relevant discussion of the reduced approximate factoriza- 
tion error that attends using DDADI in one or more space directions. 
The variety of multidimensional solution strategies derived from DDADI 
bidiagonal approximate factorization will be presented in greater detail, 
along with extension of the method to higher order, in a forthcoming 
paper 20 . 

While extensions to 3-D are not given in detail here, we note such 
extensions to a developing two level 3-D CSCM upwind scheme are pos- 
sible and will be adopted in the future. For 3-D, the inviscid terms 


T 
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similar to (6a) and (6b) will again be added for the f cross Sow direc- 
tion. For eBciency in solving each resulting marching plane, the implicit 
operators for the r, and f coordinate directions can be approsimateiy fac- 
tored using DDADI as described in reference (18) for tsro space direc 
tions. 


Two Dimensional Results 

We present results for a 45® - 15° axisymmetric transonic nozzle 
flow previously studied experimentally by Cuffd, Back and Massier 24 and 
computationally by Cline, Prozan, Serra and Shelton (all referenced in 
(24)) and ourselves 18 . In Figure 10 we show results after 10 steps of an 
early computation run at a local CFL number of 20 with the present first 
order single level scheme. Except for the addition of an error correction 
procedure 17 to counter numerical inflow boundary condition drift, a factor 
which has improved the present solution in the vicinity of the axis, the 
effectively converged results found here are the same as those given for the 
two level scheme in reference 18. (As long as the problem has a unique 
solution, the two schemes must give equivalent results since the right hand 
side difference equation sets, including boundary approximations, are the 
same.) 

For the solution given in Figure 10, we noted a very rapid rate of reduc- 
tion in residual, three orders of magnitude in ten steps. This compares 
with 60 steps given in reference (18) for the solution obtained with the 
two level scheme. The rapid convergence found in this transonic problem 
for the CSCM-S method with viscous terms provides the reasonable ex- 
pectation of similar fast results to be obtained without viscous effects. 
Thus the method in multidimensions appears to have attractive poten- 
tial for an improved transonic Euler solver as well as Navier-Stokes sol- 
ver. 
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' Nest, we present first order inviscid and viscous results for an inlet 
problem shown in Figure 11. The pressure contours for the first order inviscid 
solutions are shown in Figure 12. Figure 13 shows the first order viscous 
results. The viscous computation shows the presence of the leading edge 
shock. The flow structure compares very well with the theoretical (for the 
inviscid case) and other computational results. In Figure 14, the inviscid 
and viscous wall pressure are compared with the exact solution (inviscid). 
Figure 15 shows the convergence history of the RMS residue of all the 
consevative flow variables for the inviscid problem solved at CFL = 100 
with 4 inner iterations at every axial location. For the inviscid case, only 
forward marching was carried out and backward marching was omitted. The 
solution has converged for practical purposes at the end of the first sweep. 
The residue reaches machine accuracy in 10 iterations. In a later paper 25 , 
we will show the residual reduction versus inner iteration number in single 
sweep solutions for supersonic flow and compare results with contemporary 
PNS procedures. 


HI. Summary and Conclusions 

V 

From both theory and computational experiments -we observe that the 
new single level eigenvector split upwind method is unconditionally stable 
and practically robust for all types of flow. Thus the method removes 
the restrictions suffered by contemporary PNS methods and the resulting 
need to use other complimentary methods. As a result, we conclude, the 
present method will permit the solution of complete flowfields about complex 
configurations in a comparatively convenient, efficient and consistent manner. 
The new method, in addition to requiring only half the dependent variable 
storage, has been found in test problems to obtain converged solutions 
in fewer iterations and less net work - through improved treatment of 
nonlinear effects - than the linearized implicit methods from which it evolved. 
For flows with subsonic zones or boundary layers with strong upstream 
influence, the use of the present alternating direction implicit procedure 
that is related to symmetric Gauss-Seidel relaxation greatly accelerates 
convergence (by order the number of grid points in a subsonic zone) relative 
to global iterations (and other schemes) that repetitively sweep only in 
the predominantly streamwise direction. Hence for flows with elliptic 
influence in the flow direction, we conclude that such unidirectional sweep 
strategies, regardless of other factors, will not result in very effective solution 
algorithms. 
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Figure 1. 
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Shubin's diverging nozzle supersonic flow solution developed in 
one forward sweep from supersonic initial data. 
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Figure 2. continued 
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Figure 3. continued 
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Figure 5. Solution to subcritical nozzle problem after 15 global 
iterations with forward marching sweeps only. 
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Figure 7. Convergence history of RMS error based on the exact steady 
solution for the subcritical nozzle flow solved with forward 
sweeps only. 
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Figure 8. Convergence history comparison betwem the full 
block tri diagonal CSCM Solver and ths CSC'i-S 
method. The solutions have reached in RMS error 
less than 1 x 10“ 5 . 
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Figure 10. Mach number contours in vicinity of the throat for the axisym- 
metrlc transonic nozzle problem solved in ten global iterations 
Other computed results of Cline (dashed line) and Prozan (chain* 
dot) and experiment of Cuffel, et al (symbols) are discussed 
in references 18 and 24. 
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Figure 12. Pressure contours from the first order inviscld solution 
for the inlet problem after 10 Torward sweeps with a 
26 by 26 unstretched grid. 
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Figure 13. Pressure contours from the first order viscous solution 
for the inlet after 20 global sweeps with a 51 by 51 
stretched mesh* Note the leading edge 3hock» 
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Figure 14. Wall pressure comparison. Line - exact inviscid solution, 
square - first order inviscid solution, dot - first order 
viscous solution 
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Figure 15. Convergence history of the RMS of the residuals for the 

inviscid first order inlet problem with 4 inner iteration 
per global 3weep. Note at the end of first sweep th« 
residue ha3 reduced and the solution has converged for 
all practical purposes. 
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